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A PARALLEL METHOD FOR SOLVING POISSON 
EQUATIONS WITH DIRICHLET DATA USING LOCAL 
BOUNDARY INTEGRAL EQUATIONS AND RANDOM 
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Abstract. In this paper, we will present a new approach for solving Poisson equations 
in general 3-D domains. The approach is based on a local computation method for the 
DtN mapping of the Poisson equation by combining a deterministic (local) boundary inte- 
gral equation method and the probabilistic Feynman-Kac formula of PDE solutions. This 
hybridization produces a parallel algorithm where the bulk of the computation has no need 
for data communications. Given the Dirichlet data of the solution on a domain boundary, a 
' local boundary integral equation (BIE) will be established over the boundary of a local region 

formed by an hemisphere superimposed on the domain boundary. By using a homogeneous 
Dirichlet Green's function for the whole sphere, the resulting BIE will involve only Dirichlet 
\ data (solution value) over the hemisphere surface, but both Dirichlet and Neumann data 

■ over the patch of the domain boundary intersected by the hemisphere. Then, firstly, the so- 

lution value on the hemisphere surface is computed by the Feynman-Kac formula, which will 
be implemented by a Monte Carlo walk on spheres (WOS) algorithm. Secondly, a bound- 
ed ■ ary collocation method is applied to solve the integral equation on the aforementioned local 
patch of the domain boundary to yield the required Neumann data there. As a result, a local 
method of finding the DtN mapping is obtained, which can be used to find all the Neumann 
data on the whole domain boundary in a parallel manner. Finally, the potential solution 
in whole space can be computed by an integral representation using both the Dirichlet and 
Neumann data over the domain boundary. 

. Key words. DtN mapping, last-passage method, Monte Carlo method, WOS, Poisson 

equations 
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\ both engineering applications and mathematical study of elliptic PDEs. In 

the electrostatic potential problems, the surface charge distribution Os on the 
surface of a conductor fi, as required in the capacitance calculation of 
conductive interconnects in VLSI chips, is exactly the normal derivative of the 
electrostatic potential u as implied from the Gauss's law for the electric field 
E = — Vu, i.e., 
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1. DtN mapping and solutions of potential equations. The DtN 

mapping between the Dirichlet data (solution value) and the Neumann data 
(the normal derivative of the solution) of a Poisson equation is relevant in 
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= E-n|gf^= - — loo- (1.1) 

On the other hand, the DtN mapping also plays an important role in 
the study of the Poisson equations. As the inhomogeneity right-hand side 
of a Poisson equation is usually known, we could use a simple subtraction 
technique to reduce the Poisson equation to a Laplace equation, but with a 
modified boundary data. Therefore, in the rest of this paper we will present 
our method for the Laplace equation in a domain Q where a general Dirichlet 
data is given on the boundary dQ. If we are able to compute the Neumann 
data for the given Dirichlet data, namely the following DtN mapping: 

du 

DtN: u\an ^ Trldn, (1-2) 
on 

then, the solution u(x) at any point x in the whole space can be found simply 
by the following integral representation. 



u(x)=/ G{^,y)^^dsy- [ ^^p^u{y)dsy, x G M^f^, (1.3) 
Jdn c'ny Jq^ Criy 

where G(x, y) is the fundamental solution to the Laplace operator, namely, 

G(x,y) = ^^^. (1.4) 
Att |x — y| 

A similar NtD mapping from Neumann data to Dirichlet data can also be 
defined if the Neumann data yields a unique solution to the PDE. In either 
case, with both Dirichlet and Neumann data at hand, the solution of a Laplace 
equation can be obtained by the representation formula in (jl.3p . 

In addition to the electrostatic potential problem in the capacitance calcu- 
lation, the solution of the Poisson equation or Helmholtz equations is also the 
main component of projection-type methods for incompressible flows [5]|27j. 
which usually involve the solution of a Helmholtz equation for the velocity field 
and a Poisson equation for the pressure field. Therefore, by finding the DtN 
or NtD mapping of the relevant elliptic PDE solutions in an efficient manner, 
we could produce fast numerical methods for many electrical engineering and 
fluid mechanics applications. 

For the capacitance problems of conductors, boundary element method 
(BEM) or finite element method (FEM) have been used as electrostatic field 
solvers to compute the charge density, for example, the indirect BEM Fast- 
Cap [2i][25], the direct BEM QMM-BEM [29], hierarchical extractors Hi- 
Cap and PhiCap [26] [28], and the parallel adaptive finite element method the 
ParAFEMCap [7], etc. BEMs [2] need to discretize entire conductor surfaces, 
sometimes even the dielectric interfaces, into small panels, and construct a lin- 
ear system by method of moments or collocation methods to be solved. These 

2 



deterministic methods are highly accurate and versatile, but are global, i.e., 
even if the charge density at only one point is required, a full linear system has 
to be constructed and solved. In general, the resulting linear algebraic sys- 
tems are solved by iterative methods such as the multi-grid methods [3] or the 
domain decomposition methods |30j . either as a solver or pre-conditioner. For 
integral equation discretization, the fast multipole method (FMM) [13] can be 
used in conjunction with a Krylov subspace iterative solver. All these solvers 
are 0{N) in principle and iterative in nature, and require expensive surface 
or volume meshes. The parallel scalability of these solvers on large number 
of processors poses many challenges and is the subject of intensive research. 
In comparison, the method proposed in this paper is intrinsically parallel and 
has the potential of high parallel scalability. 

Due to the limited memory in computers and computational time, some- 
times it is impractical to obtain solutions by global methods for many engineer- 
ing problems such as modern VLSI chips with millions of circuit elements. In 
contrast, random methods can give local solutions, which have been applied in 
the area of chip design industrials. For instance, the Rapid 3D and QuickCap, 
as the chip industry's gold standards produced by the leading EDA companies 
Synopsys and Magma respectively, are both random methods. The key ad- 
vantage of the random methods is their localization. For example, QuickCap 
[8] [9] can calculate the potential or charge density at only a point locally with- 
out finding the solution elsewhere. Usually, random methods are based on the 
Feynman-Kac formula and the potential (or charge density) is expressed as a 
weighted average of the boundary values [T7]. 

The Feynaman-Kac formula |12] relates Ito diffusion paths to the so- 
lution u(x) of the following general elliptic problem 

L{u) = Yhi{^)^^ + W^^^ = ^^^^ 

i=l * iii = l * 

u\qq = 0(x), X € dQ. (1.5) 

If Xt{u}) is an Ito diffusion defined by the following stochastic differential 
equation 

dXt = h{Xt)dt + a{Xt)dBt, (1.6) 

where Bt{oj) is the Brownian motion, [ajj] = ^a(x)a^(x), then, the following 
Feynman-Kac formula gives a probabilistic solution for (jl.Sp as 

u{^) = E'-{cPiXrJ + E^[ f{Xt)dt], (1.7) 

Jo 

where the expectation is taken over all sample paths Xj=o(f^) = x and tq is 
the first hit time (or exit time) of the domain 17. This representation holds for 
general linear elliptic PDEs. However, for the purpose of this paper, we will 
only consider the Laplace equations. 
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For the Laplace equation (/ = 0), the Ito diffusion is just the Brownian 
motion and the solution can be simply rewritten in terms of a harmonic mea- 
sure, which measures the probability of the Brownian paths hitting a given 
area on the boundary surface, 

u(x) = E-(</)(X.J) = / Hy)dt^f„ (1.8) 
Jan 

where 

= P^uj\Xrn{uj) e F,Xo(w) = x}, F C dn, x e f). (1.9) 

The harmonic measure can be shown to be related to the Green's function 
^(x, y) of the Laplace equation in the domain Q with a homogeneous boundary 
condition, i.e., 

-A5(x,y) = 5(x - y), x e 1), 
<?(x,y)|xean = 0. (1.10) 

By comparing (ll.Sp with the following integral representation of the solu- 
tion of the Laplace equation with the Green's function g(x,y), 

f \ /■ ,/ ^ (9g(x,y) 

= - / <p(y)^ dsy, 

Jan Cny 

we conclude that the hitting probability, now denoted as p(x, y)ds. 
dsy]), has the following connection to the Green's function |10j . 

p{x,y) = . 1.12 

dUy 

Therefore, if the domain is a ball centered at x where a path starts, we 
have a uniform probability for the path to hit the surface of the ball. This 
fact will be a key factor in the design of random walk on spheres (WOS), 
which allows us to describe the Brownian motion and its exit location without 
explicitly finding its trajectory. Instead, a sequence of walks or jumps over 
spheres will allow the Brownian path hitting the boundary (for practical 
purpose, within an absorption e-shell as proposed in \2T\). Specifically, as 
indicated by (jl.l2p . the probability of a Brownian path hitting on a spherical 
surface is given exactly by the normal derivative of the Green's function of the 
sphere (with homogeneous boundary condition). Therefore, if we draw a ball 
centered at the starting point x of a Brownian path, it will hit the spherical 
surface with a uniform probability as long as the ball does not intersect with 
the domain boundary dQ. So, we can make a jump for the Brownian particle 
to xi, sampled with a uniform distribution on the spherical surface. Next, a 
second ball now centered at xi will be drawn, not intersecting with the domain 
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(1.11) 




Fig. 1.1. WOS sampling of Brownian pathes 

boundary dfl, the Brownian particle can make a second jump to X2 on the 
surface of the second ball. This procedure (as illustrated in Fig ll.H termed 
as WOS) |22j|14j[25]. is repeated until it hits the boundary of Q (within the 
e-shell of absorption) where it is denoted as x^-j^ and the value of the boundary 
data (p{xr^) will be recorded and eventually all such data will be used to 
compute the expectation in (II. Sp . 

In real applications, due to the relation between the Green's function 
g{x, y) of a domain and the hitting probability, Green's Function First Passage 
(GFFP) methods for shapes other than spheres such as rectangles in softwares 
including QuickCap [8] [9] have been used to find capacitances of conductors 
in interconnect layouts, which are generally of rectangular shapes. 

The Feynman-Kac formula allows a local solution of the PDE, and fast 
sampling techniques of the diffusion pathes with the WOS methods are avail- 
able for simple PDEs such as Laplace or modified Helmholtz equations. How- 
ever, it is impractical to use the probabilistic formula to find the solution of 
these PDEs in whole space as too much sampling will be needed. In this 
paper, we will propose a hybrid method for computing the DtN mapping by 
combining the probabilistic Feynman-Kac formula and a deterministic local 
integral equation over domain boundary dft. The hybrid method will allow us 
to get the Neumann data efficiently over a local patch of the domain boundary, 
which will result in a simple parallel method for solving the complete potential 
problems in general 3-D domains. 

We will present the work in the following stages. Firstly, we will review 
a last passage method in [H] which calculates the Neumann data (charge 
distribution) at one single point over a flat surface when the Dirichlet data is 
a constant. Even though this is a very limited case for the DtN problem, it 
demonstrates some key issues and difficulties in how to use the Feynman-Kac 
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formula and the WOS in finding the Neumann data. Secondly, we will present 
our hybrid method, which allows the calculation of Neumann data for a general 
Dirichlet data on the flat surface. Thirdly, we extend the hybrid method 
to calculate the Neumann data over a patch of the boundary for arbitrary 
Dirichlet data and curved boundary. Numerical tests will be presented to 
show the accuracy and potential of the proposed methods. Finally, conclusion 
and discussions for open research issues and parallel aspect of the proposed 
method will be given. 

2. Finding the Neumann data at one point over a flat boundary. 

2.1. Review of a last-passage algorithm for charge density. In this 
subsection, we will review the last-passage Monte Carlo algorithm proposed 
in [13] for charge density, namely the Neumann data, at one point on a flat 
conducting surface. 

For a flat portion of the boundary d^l of a domain Q = {z < 0} in the 3-D 
space held at a constant potential, we like to compute the charge density at a 
point X € do,. In the last-passage method, a hemisphere is constructed with a 
radius a centered at x as shown in Fig. 12.11 The upper hemispherical surface 
outside Q is denoted as F and the 2-D disk of radius a centered at x from the 
intersection of the hemisphere and the conducting boundary dO, is denoted as 

Sa = 5a(x). (2.1) 

From the discussion in Section 1 on the equivalence between the electro- 
static potential and diffusion problems, the quantity u(x) = 1 — u{x + e) can 
be interpreted as the average value of u = on disk Sa, i.e., 

v{x) = 0, xe5a (2.2) 

and V = 1 at infinite (or on a infinitely large sphere). Therefore, we have the 
following probabilistic expression for v{x-\-e), also viewed as the probability of 
a Brownian particle near the conducting surface d^} starting at x + e diffusing 
to infinity without ever coming back to the conducting surface |14] |22j : 

v{x + e) = 1 - u{x + e) = - j g{x + e, y)pyoodsy, (2.3) 

where Pyoo is the probability of a Brownian particle starting at y and diffusing 
to infinity without ever coming back to the conducting surface dO,, thus, py^o = 
if y € 5*0. In (12. 3p . the integral over F expresses the Markov property of 
the diffusing particles from x + e to the infinity with an intermediate stop 
on F. Specifically, ^(x + e,y) gives the probability of a Brownian particle 
starting from x + e hitting the boundary of F, which is given by (jl.l2p via a 
homogeneous Green's function for the hemisphere over Sa, namely, 

da 

?(x + e,y) = -^(x + e,y), (2.4) 

oily 
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Fig. 2.1. Last passage for finding the Neuamann data at one point 

and 5(x + e, y) is defined in (jl.lOp for tlie hemispliere, whose analytical form 
can be obtained by an image method with respect to spherical surface first, 
then to the half plane z = 0, resulting in the use of 3 images. Specifically, we 
have 

5(x,xs = — ^ \ + ir'\ r + :ri i+:ri T' (^-^^ 

47r|x — Xs| 47r|x — Xk| 47r |x — Xs| 47r |x — X]^| 

where in spherical coordinates the source location is Xg = {ps,0s,4's), the Kelvin 
image location with respect to the sphere is x^ = (a^/ps, ^g, (/)s),and their 
mirror image locations with respect to the plane z = are Xg = (/9s,vr — 
6*8, <?!)s), xj- = {a? / ps,'K — 9s, (ps), respectively. Meanwhile, the corresponding 
charges are gk = —a/ps, Qs = —1, and = a/ps, respectively. 

Now, to get the charge distribution as (normal derivative), we use the 
relation in (jl.ip . we have 

as = -limn^+e • E(x + e) = lim^^^^^i^ = (2.6) 



and 

c^'u(x) f 

— = h(x + £,y)pyoodsy = T,i^p, (2.7) 

where a shorthand Slp is introduced for the integral over T for latter use, and 



^2 

The weight function /i(x, y) can be analytically computed for the hemisphere 



3 cos 6* 

/ix>y) = 7^ 3-) (2-9) 
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where is the angle between the two normal vectors n'^ and on the boundary 
r as shown in Fig. 12. li 

Next, we only need to compute pyoo which is the probability of a Brownian 
particle starting from y G T and diffuse to infinity without ever returning to 
the conductor surface dn. Due to the homogeneity of the Brownian motion 
in the external domain = {z > 0}, the WOS in Section 1 can be used to 
calculate this probability. The integral in ()2.7p over T could be approximated 
by a Gauss quadrature as both /i(x, y) and pyoo can be considered as smooth 
functions of y € F. Nonetheless, in [14], the integral Slp is computed by first 
distributing A'^ particles at locations over F based on a distribution density 
derived from (j2.9p . and then starting a Brownian diffusion path from each of 
those locations. The number of paths which will diffuse to infinity (in practice, 
when it hits a very large ball) is recorded as A^inf , then, we have the following 
estimate 



The key equation in the last-passage algorithm is (j2.7p . which is based on 
(j2.3p provided that the potential solution v{x.) = 0,x € 5a on the conductor 
surface as indicated in (12. 2p . Therefore, for general non-constant Dirichlet 
boundary data, the last passage method will not be applicable. In fact, the 
charge density at x will be influenced by potential value on all domain bound- 
ary. 

2.2. BIE-WOS Method: Combining a BIE and the Monte Carlo 
WOS method. For the last-passage method discussed above, the algorithm 
(j2.7p is obtained by the isomorphism between the electrostatic potential and 
diffusion problems. The limitation of the last passage method is that it is only 
applicable to the situation of constant Dirichlet data and a flat boundary. In 
this section, we will adopt a different approach based on a boundary integral 
equation (BIE) representation of the charge density (the Neumann data) on 
the surface at a given point using potential over a small hemisphere, the latter 
can be then computed by the random WOS method. As a result, this new 
approach, a hybrid method of deterministic and random approaches, will be 
able to handle general variable Dirichlet boundary data, and later in Section 
3 also be extended to curved boundaries. 

Let us denote by il^ the domain formed by the hemisphere of radius a 
centered at x over the flat boundary boundary Sa as in Fig. 12.11 by applying 
the integral representation (jl.Sp of Laplace equation with the afore-mentioned 
Green's function ^(x, y ) in (12. Sp for the domain with a homogeneous Dirich- 
let boundary condition. Due to the zero boundary value of the Green's function 
5(x,y), we have 
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(2.11) 



where F again is the surface of the upper hemisphere and Sa is the disk of 
radius a centered at x. In order to obtain the normal derivative of u at x, 
we simply take the derivative with respective to x' along the direction nx as 
x' approaches x and obtain the following representations involving a hyper- 
singular kernel, 

-^n(x) = - hm / ^p^^u{y)ds, x G 5,. (2.12) 

The integral expression for ^^u(x) involves two integrals, one regular 
integration over the upper hemisphere T denoted as 



X'— >-x 



d 



where (12. 9p has been used in the second equality, and another hyper-singular 
integral over the disk Sa denoted as 

g^g(x',y 
J duydn 

and we have 

-^n(x) = Si + S2. (2.15) 

Equation ()2.15p will be the starting point for the proposed hybrid method. 
In computing the integral Si, say by a Gauss quadrature over the hemisphere 
surface, we will need the potential solution u{y) for y € F and this solution 
will be readily computed with the Feynman-Kac formula (jl.Sp with the WOS 
as the sampling technique for the Brownian paths. On the other hand, the 
singular integral S2, with appropriate treatment of the hyper-singularities to 
be described in detail in the numerical test section, can be calculated directly 
with the given Dirichlet boundary data w(y),y G Sa- Therefore, an algorithm 
using (j2.12p involves the hybridization of a random walk on spheres (WOS) 
and a deterministic boundary integral equation (BIE), which is termed the 
BIE-WOS method. 



Remark: In comparing the last-passage method ()2.7p and the BIE-WOS 
method (j2.15p . the former uses the relation between the Brownian motion of 
diffusive particles and electric potential from charges on a conducting surface 
to arrive at an expression for the surface charge density based on ()2.7p . The 
BIE-WOS uses an hyper-singular boundary integral equation to get a similar 
expression in ()2.15p , which has an additional contribution from the variable po- 
tential on the charged surfaces (the integral term S2). Both methods use WOS 
for particles starting on the hemisphere, however, at different locations. The 
last-passage method proposed in (|14j) initiates particles walk starting from 
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Fig. 3.1. Setup of the BIE- WOS method for finding Neumann data on a patch S C dfl. 

positions all over the hemisphere sampled using a probability given by (j2.9p 
while BIE- WOS initiates many particle walks starting from selected Gauss 
quadrature points (up to 30 x 30 in our test problems). Numerical results will 
show that for problems suitable for both methods, the total number of particle 
walk paths and the accuracy and computational costs are comparable (refer 
to Test 4 in Section 4.1.3). 

3. Finding Neumann data over a patch of general boundary. In 

this section, we will extend the BIE- WOS of last section to the case of general 
Dirichlet boundary data and curved domain boundary. To achieve this goal, 
we will superimpose an hemisphere over any selected portion of the boundary 
dil. and denote the intersection portion of the domain boundary by S and 
the surface of the hemisphere outside the domain 0, still by T and the region 
bounded by S and T is denoted as Qs (see Fig. 13. ip . Now let G(x, y) be the 
Green's function of a whole sphere with an homogeneous boundary condition 
and G(x, y) can be easily obtained by one Kelvin image charge as discussed 
before. Then, the integral representation ()1.3p can be applied to the boundary 
of the domain to yield the following identity 



'u(x) 



gg(x,y) 
5nv 



u{y)dsy+ 



gg(x,y) 
5nv 



n(y) + G(x,y) 



du{y) 



d 



ds 



X e 1^5 



(3.1) 

It should be noted that the integral over F only involves the normal derivative 
of the Green's function as G vanishes on F by construction. As a result, 
only solution u{y) is needed on F while both u{y) and the normal derivative 
appear in the integral over S. As before, the solution u{y) over F will 



be computed with the Feynman-Kac formula (|1.8p with WOS and then the 
Neumann data over S can be solved from the following integral equation. 



K 



du 
dn 



(x) = 6(x), XG5, 



(3.2) 
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where 



K 



and 



6(x) 



u(x) 



+ p.v. 



du 
dn 



aG(x,y) 



G[x,y)— (is 



d. 



u{y)dsy) + 



d 



(3.3) 



u{y)dsy. (3.4) 



Integral equation (|3.2p is of the first kind which is ill-conditioned and 
may cause numerical difficulties especially when the algebraic system from 
discretization becomes large. When that happens, a well-conditioned second 
kind of integral equation can be obtained by taking normal derivative of (j3.ip . 
resulting in the following identity for x G 



d 



d: 



u(x) 



d^GiK,y) 
r 3ny3nx 



u{y)dsy+ 



a2G(x,y) 9G(x,y)an(y) 

-u{y) + 



d: 



n 



y J 
(3.5) 

Let X approaching the boundary S, we obtain the following second kind inte- 
gral equation 



1 f)ii 

(-/-!))[— ](x) = 6(x), xe5, 

2 on 



where the integral operator of a double layer potential 

D[—]{^) = f ^G'(x,y) dujy) ^^ 
dn Js 9nx dny 



(3.6) 



(3.7) 



and 



6(x) ^ - f ^^^^'l;^\ {y)dsy - p.f. / ^^^^Y^ u{y)dsy, x G 5, (3.8) 

and p.f. denotes the Hadamard finite part limit for the hyper-singular integral, 
which can be handled by a regularization technique. 

BIE-WOS Algorithm: The BIE-WOS method for the Neumann data 
over a patch S will consist of two steps: 



• Step 1: Apply the Feynman-Kac formula (jl.Sp with the WOS sample 
technique to compute the potential solution at Gauss points yjj G F, 
u{yij). Compute the right hand side function 6(x) in (13. 4p or (13. 8p by 
some Gauss quadratures. 

• Step 2: Solve the BIE ()3.2p or (j3.6p with a collocation method for the 
Neumann data ^ over S. 
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4. Numerical Results. In this section, we will present a series of numer- 
ical tests to demonstrate the accuracy and efficiency of the proposed BIE-WOS 
method for finding the Neumann data at a single point on a flat boundary and 
on a patch over a curved boundary. 

4.1. Finding Neumann data at one point on a flat boundary. 

4.1.1. Regular izat ion of hyper-singular integrals. First, let us present 
a regularization method using simple solution of the Laplace equation [16j to 
compute the hyper-singular integral in ()2.14p and ()3.8p . First, with some sim- 
ple calculations, the term T,2 of (|2.14p is found to be a Hadamard finite part 
limit of the following hyper-singular integral: 

T,2 = - lim / ^ '^\ {y)dsy = -p.f. / ^ - u{y)dsy, (4.1) 

where /? = |x — y|, x, y G Sa- The finite part (p.f.) limit of Hadamard type 
is defined by removing a divergent part in the process of defining a principal 
value (i.e. by removing a small patch of size e centered at x and then let 
£ approaching zero). For the Laplace equation considered here, we can reg- 
ularize this hyper-singularity by invoking an integral identity for the special 
solution u = (/>(x), X is fixed, namely, the integral identity (I2.12p applied to 
this constant solution results in 



= - / f^0(x)^. - lim / ^^^(x)d., X , S. (4.2) 

Subtracting (|4.2p from (j2.15p . we have a modified formula for the Neumann 
data as 

-^n(x) = S; + S'2, X € 5, (4.3) 
orix 

where S'^ and S2 are now regularized versions of Si and S2 in (j2.13p and 
(|2.14p . respectively, i.e., 

and 

where x' = x+(0, 0,e), r = sj + e^^p = [x — y[,x,y G Sa- Moreover, the 
boundary condition u(y) = </>(y), y € 5a has been invoked in (j4.5p . 
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Compared with (j4.ip . the singularity in the integral in (|4.5p has been 
weakened by the factor {4>{y) — (/"(x)) , which vanishes at x, and will be 
evaluated by a Gauss quadrature. Let us only consider the integral involving 
the singular term ^ in ()4.5p . which is denoted by II2, i.e., 

= hm / 1 (<A(y) - </.(x)) ds,. (4.6) 

Consider a circular patch A5 of radius 5 centered at x, and then can 
be split further into two integrals as follows 



1 f 1 



s 

= -^l -,{m-m)dsy + ^. (4.7) 

To estimate the term A, we apply a Taylor expansion of the boundary 
data (piy) at x 



x) = V,/.(x).p + 0(p2), (4.8) 



then, we obtain 



V(/.(x)- ^. f p ^ If 0{p 



A = lim / -^dsy + — — 5 — ds 



27r x'^xjj^ ^ 2tt 



3 ""f 



v^(x). 



2^ x'^x7o io (p2 + e')'/' 
(p2 +£2)3/2' 

<5 rjC^s^ 



+ 1™ / T^T— ^^T^/'^P 

x'— )-x 



X'— >x 



Now for all positive e > 0, we have 

as a result, the following estimate of the term A holds 

A = 0(5). (4.11) 

Finally, the regularized integral S2 will be approximated by the integral over 
Sa\As with an accuracy of 0{5) and a Gauss quadrature formula over the ring 
shaped region Sa\As: 

^*2 = -^f 4 - + 0{S). (4.12) 
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4.1.2. Gauss integrals over the hemisphere F and Sa\As and WOS. 
To compute the integral T,[, we use NgixNgi Gauss points over the hemisphere 
surface T 

S; -J2u^,u^j^(a'sme,)^ (^(y,,,) - 0(x)) , (4.13) 

where 

Qi = \{ii^ l),^i = + l),yM = (a,0i,(^j)> (4.14) 

and LOi and ■^j, 1 < ^ < A'gi are the Gauss quadrature weights and locations, 
respectively, ^(a^sin^j) is the surface element in the spherical coordinates. 

Now, each of the solution value ■u(yjj), yjj- S F will be obtained by the 
Feynman-Kac formula (jl.Sp with Npath Brownian particles all starting from 
yi,j, namely 

<y^,,)^^r- E 't'iek), (4.15) 

^^P^^^path k=l 

where e^ is the location on dVt where a path terminates. 

The total number A^path-bic-wos of Brownian particles needed in the BIE- 
WOS method will be 

^path — bic— wos 

X Ngi X Npath- (4.16) 

Next, the integral S2 in (j4.12p will be computed with another Ng2 x Ng2 
Gauss quadrature over the ring shaped region Sa\-^-& with an error of 0{5) in 
addition to the error from the Gauss quadrature. 

4.1.3. Numerical tests. In this section, we will present several numeri- 
cal tests to demonstrate the accuracy and efficiency of the proposed BIE-WOS 
for finding the Neumann data at a given point over a flat boundary for general 
Dirichlet boundary data. For comparison, we also implement the last-passage 
Monte Carlo method in [14]. For accuracy comparison, the charge density is 
calculated with the FastCap, an open-source code developed in MIT [2l] for 
3-D capacitance extraction tool in industry and academia. The Fastcap is an 
indirect BEM, accelerated by the fast multipole method (EMM), and its lin- 
ear system is solved by a conjugate gradient method. For the case of complex 
potentials on the surfaces, we also implemented a direct BEM (DBEM) |29j . 

• Test 1- Charge densities on a planar interface between two dielectric half 
spaces 

As shown in Fig. 14. H the whole space is divided by a planar interface 
between the two dielectric domains, and the dielectric constants are eo and 
ei in the upper and lower domain, respectively. A charge q is located at 
r^ = (0, 0, —h) and then the potential in the upper space is given by 
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Fig. 4.1. Potential above a half-space 



Table 4.1 

CHARGE DENSITY ON THE HALF PLANE INTERFACE WITH DIFFERENT RA- 
DIUS 



a 


Last-passage 


Bie-WOS 


analytical 
solution 




err% 






s; + 


err% 


0.1 


0.698543 


-2.38 


0.69884 


0.018777 


0.717612 


0.29 


0.71554 


0.2 


0.677996 


-5.25 


0.67784 


0.037515 


0.715355 


-0.03 


0.5 


0.622949 


-12.94 


0.62146 


0.093054 


0.714517 


-0.14 


0.7 


0.586721 


-18.00 


0.58432 


0.128971 


0.713287 


-0.32 


1.0 


0.534695 


-25.27 


0.53659 


0.179973 


0.716562 


0.14 



u{r) = 1 r, q = ■ q, 4.17 

47reo|r-r^| gq + ei 

and u{r) satisfies the Laplace equation V^ii(r) = 0, z > with a variable 
Dirichlet data on the boundary z = 0. 

The charge density at the point x = (0.5,0,0) by the last-passage method 
and BIE-WOS with various radius a of the hemisphere are listed in Table 
14.11 Li the last-passage method, the total number of the sampling paths 
AT = 4 X 10^ In the BIE-WOS, the number of Gauss points Ngi x Ngi = 20 x 20 
for the hemisphere and Ng2 x Ng2 = 20 x 20 for the integral on the 2-D disk 
Sa, and starting from each Gauss point on the hemisphere, the number of 
the sampling paths is Npath = 10^. Therefore, the total number of paths for 
the BIE-WOS method is also 4 x 10^. In both methods, the thickness of the 
absorption layer e for the WOS method is taken to be 10~^, and the step 
threshold for going to infinity Ngtep = 300. 

From Table 14. H we can see that when the radius becomes larger, the 
relative errors of the last-passage method grows and even up to —25.27%. It 
shows that when the potential Dirichlet data on the disk Sa is not constant, 
the last-passage method is not applicable. The variable potential inside the 
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Table 4.2 
ACCURACY OF WITH AND 



6/a 


S2 calculated by Ng2 x Ng2 Gauss Quadrature 


4x4 


err% 


6x6 


err% 


10 X 10 


err% 


20 X 20 


err% 


10-1 


0.09659 


3.800 


0.08983 


-3.462 


0.08949 


-3.833 


0.08949 


-3.8351 


10-^ 


0.10042 


7.916 


0.09306 


0.001 


0.09273 


-0.347 


0.09273 


-0.3492 


10-^ 


0.10083 


8.352 


0.09335 


0.314 


0.09302 


-0.033 


0.09302 


-0.0346 


10-4 


0.10087 


8.397 


0.09337 


0.345 


0.09305 


-0.002 


0.09305 


-0.0035 


10"^ 


0.10087 


8.402 


0.09338 


0.348 


0.09306 


0.001 


0.09305 


-0.0003 


10-^ 


0.10087 


8.402 


0.09338 


0.348 


0.09306 


0.001 


0.09305 





disk Sa at the bottom of the hemisphere will influence the charge density at x 
to be calculated. In contrast, the BIE-WOS includes such influences as shown 
in the results, and most importantly, it is independent of the radius a, for its 
maximal relative errors is less than 0.32% when the radius ranges from 0.1 to 
1.0. 

Table 14.21 lists the accuracy of the de-singularized S2 in (j4.12p with differ- 
ent values 6 and number of Gauss points Ng2 x Ng2, where the location of the 
sought-after density is at (0.5,0,0). The result of Ng2 x Ng2 = 20 x 20 with 
6/a = is taken as the reference value for Sg. Table H?2] shows the conver- 
gence speed of S2 as 6/a goes to zero and number of Gauss points increases. 
It can be seen that when the number of the Gauss point is large enough, for 
example 20x20, the relative error is on the order ol d/a, verifying the estimate 
in (fri2]) . 

• Test 2: Four rectangle plates with a piecewise constant potential distribution 

A 3-D structure with four rectangle plates is depicted in Fig. 14.21 where 
the length, width and thickness of all four unit plates are Im x Im x 0.01m. 
First, we set the potential of plate II to IV and other three (I, III and IV) 
to OV, and compute the charge density at the point ^(—0.2273,0.2273). The 
results of all four methods are listed in Table 14.31 taking the results by the 
FastCap as the reference where each side of the plates is discretized into 99x99 
panels. The DBEM uses a discretization with 11x11 panels on each side, its 
relative error is -0.46%. 

Both last-passage method and BIE-WOS run with various radius a of the 
hemisphere, and the parameters are the same as in Test 1. In this case, the 
integral is related to the area of the intersecting area between the disk Sa 
and the plate I, III and IV, we just compute it directly by the quad function 
in Matlab, instead of Gauss quadratures. 

Note that the potential on the boundary here is piecewise constant. 
Therefore, in the last-passage method, charge density should be computed, 
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Fig. 4.2. Four plates at different potentials 



Table 4.3 

CHARGE DENSITY OF A FOUR UNIT PLATES STRUCTURE WITH DIFFERENT 
RADIUS 



a 


Last passage 


BIE-WOS 


DBEM 




err% 






s; + S2 


err% 


value 


err% 


0.1 


2.6084 


0.05 


2.6051 





2.6051 


-0.07 


2.595 


-0.46 


0.2 


2.6026 


-0.17 


2.6051 





2.6051 


-0.07 


0.2273 


2.6099 


0.11 


2.6064 





2.6064 


-0.02 


0.3 


2.5252 


-3.14 


2.5178 


0.0892 


2.6070 


-0.00 


Fastcap 


0.5 


1.9698 


-24.44 


1.9692 


0.6330 


2.6022 


-0.19 


2.607 


0.7 


1.5779 


-39.48 


1.5784 


1.0271 


2.6055 


-0.06 



instead of by (I2.10p . by the following formula: 

3 iVinf + Nj + Niii + Njy 
Slp = , (4.18) 

jVpath-LP 

where A'^inf , Ni, Nm, and Niv represent the number of particles which finally 
go to infinite, plate I, III, and IV, respectively. A^path-LP denotes the total 
number of Brownian pathes starting from the hemisphere T. 

From Table 14.31 we can see that when the radius a < 0.2773, the disk 
Sa is totally inside the plate II, and the last-passage method is correct with a 
maximal relative error less 0.17%. However, once Sa becomes larger and covers 
areas of plates with different potentials, the relative errors of the last passage 
method increases and even up to —39.48%. In comparison, the BIE-WOS 
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Fig. 4.3. Convergence of BIE-WOS vs number of Brownian pathes and Gauss points 

maintains its accuracy insensitive to the radius a with a maximal relative errors 
less than —0.19% as the radius varies from 0.1 to 0.7. This again confirms the 
fact that the last-passage method of [H] is designed for conducting surfaces 
(i.e., constant potential), not for surface of variable potentials. Therefore, it 
should not be used when the disk Sa includes regions of different potentials. 

In conclusion, for a general variable potential, the last passage method is 
limited while the BIE-WOS does not suffer from the constrain of a constant 
boundary potential. 

• Test 3: Four rectangle plates with a complex potential distribution 

To further emphasize the point raised above in Test 2, we set the four 
plates with a complex potential distribution as: 

(/)(x,y) = sinmxsinny. (4.19) 

To obtain an accurate result, the last-passage method will require increasingly 
smaller radius a for ever larger m and n to achieve an (approximately) constant 
potential within the disk Sa- 

The charge density at the point (-0.5, 0.5) by the last-passage, BIE-WOS 
and DBEM are shown in Table WM We take the result of DBEM with 17x17 
panels on each plate as the reference solution. All other parameters in BIE- 
WOS and last passage methods are same as in the previous case. From Table 
14.41 we can see that the BIE-WOS is more accurate. 

The relative errors versus the number of Gauss points and the WOS paths 
are shown in Fig. The BIE-WOS result of iVgiX Ngi = 20 x 20, Ng2X 
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Table 4.4 

CHARGE DENSITY OF A FOUR UNIT PLATES STRUCTURE WITH COMPLEX 
VOLTAGES IN DIFFERENT RADIUS 



a 


Last passage 


BIE-WOS 


DBEM 




err% 








err% 


value 


0.1 


-0.4522 


3.92 


-0.4454 


-0.008617 


-0.4540 


3.54 


-0.4707 


0.2 


-0.4442 


5.62 


-0.4444 


-0.01722 


-0.4616 


1.93 


0.3 


-0.4369 


7.17 


-0.4362 


-0.02579 


-0.4620 


1.85 


0.4 


-0.4288 


8.90 


-0.4278 


-0.03433 


-0.4621 


1.82 


0.5 


-0.4203 


10.7 


-0.4202 


-0.04280 


-0.4630 


1.62 



Ng2 = 10 X 10 and Npath = 2 x 10^ and a = 0.5 is taken as the reference 
solution. From Fig. 14.31 we can see that when the number of the Brownian 
paths Npath is larger than 10^ at each Gauss point, the BIE-WOS result with 
10 X 10 Gauss points will reach an accuracy about 1% in relative error. 

• Test 4-' CPU time comparison with last-passage method 

For both the last-passage and BIE-WOS methods, the CPU time is ex- 
pected to be linear in terms of the total number of random paths. We demon- 
strate this fact with of a thin circular disk with radius b case in 3-D 
space [13] as shown in Fig. 14. 4[ From the analytic result of the charge 
density on the disk is: 

For a given relative error tolerance on the charge density at (-0.5,0,0) , the 
CPU time comparison of both methods versus the number of random paths 
are listed in Table H31 We take the radius a = 0.4 for 6 = 1 for the radius 
of the thin disk, and the analytic charge density is it(0.5) = 0.735105. From 
Table 14.51 we can see that the CPU times are indeed in proportion to total 
number of random paths for both methods for a comparable accuracy. Though 
the integral S2 of the BIE-WOS method in this case is obviously zero, we still 
evaluate it just as for a general variable potential and the CPU time of S2 is 
included in the CPU time of the BIE-WOS method in Table 14.51 It is noted 
that the CPU time in computing the integral for all cases are insignificant 
at about 0.012 second for a 20x20 Gauss quadrature. 

4.2. Finding Neumann data over a patch of a curved boundary. 

Next, to test the BIE-WOS method for a curved boundary, we compute the 
DtN mapping on a big sphere as shown in Fig. 14.41 with a radius R = 3. A 
point charge g = 1 is located at the central point O and the analytic result 
for the potential is then known. To compute the Neumann data over a local 
patch S around the point o = (0, 0, 3) on the big spherical surface, a small 
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r b 

Fig. 4.4. Finding the charge distribution over a disk m 3-D 
Table 4.5 

THE RELATIVE ERRORS AND THE CPU TIMES COMPARISON ACCORDING 
TO THE NUMBER N OF RANDOM PATHS 



Last passage (LP) 


BIE-WOS 


^path-LP 


Slp 


err% 


cpu 
time(s) 


-'^path— bic— wos 




err% 


cpu 
time(s) 




0.69975 


-4.81 


32 


10^ • 100 = 10^ 


0.68888 


-6.29 


30 




0.73253 


-0.35 


331 


10^ • 1000 = 10^ 


0.73960 


0.61 


307 


4-10^ 


0.73743 


0.32 


1325 


20^ • 1000 = 4-10^ 


0.73441 


-0.09 


1218 



sphere with a radius a = 1 is superimposed over the point o. The local patch 
S is discretized with a triangular mesh as shown in Fig. 14.51 

The BIE equation of (j3.2p is solved by a collocation boundary element 
method. When the collocation point is not inside an integration panel of the 
BEM, a simple Gauss quadrature method is used. For collocation point inside 
an integration panel, both weak and strong singularities will occur, however, 
they can be regularized by a local polar transformation technique and a 20 x 20 
Gauss quadrature will then be used. For the integrals on F, a 30 x 30 Gauss 
quadrature will be used. The potential n(y) on F is first computed, by the 
Feynman-Kac formula and the WOS method with lO'^ Brownian paths, on a 
regular grid, which is generated by evenly discretizing the spherical surface 
along the polar and azimuthal angles. The value u(y) at other locations on F 
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Fig. 4.5. Left: BIE-WOS setting for finding Neumann data over a patch of a bigger 
sphere; right: the mesh over the patch 



0.2 0.3 



X XX 

X jc ^ >^ ^ 



Fig. 4.6. Accuracy of the Neuamann data by the BIE-WOS over the patch S,r < 0.7a, 
a is the radius ofV. 



as required by numerical quadratures will be interpolated using the values on 
the grid points. 

The relative errors at the center of triangular panels on S are shown in 
Fig. I4.6| where x-axis means the distance between the triangle center to point 
o. From Fig. 14.61 we can see that for the panels close to point o, i.e. r < 0.7a, 
the maximal relative error is less than 1.25%, which will be accurate enough 
for most engineering applications. It is noted that due to the sharp corner 
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edge singularity of the domain where the hemisphere and d^l intersect, 
the collocation BEM will lose some of its accuracy, which limits the region 
where acceptable accuracy of the BEM solution can be used for the sought- 
after Neumann data. This well known problem in singular boundary elements 
usually is addressed with graded mesh near the edge singularity [6] [l] fl9j and 
still an active research topic in BEM methods [1]. A resolution of this edge 
singularity can increase of the region of useful BEM solution in the BIE-WOS 
algorithm and can be incorporated into the algorithm. As discussed in the 
last section, as the boundary d^l will be covered with an overlapping patches 
Si, the loss of the accuracy of the BIE solution near the edge of each patch 
will not hinder the use of the BIE-WOS method. However, any improvement 
of the BEM near the edge will reduce the total number of patches to cover the 
boundary, thus reducing the total cost. 

5. Conclusion and discussions. In this paper we have proposed a local 
BIE-WOS method which combines a local deterministic singular boundary 
integral equation method and the Monte Carlo WOS algorithm to find the 
Neumann data on general surfaces given Dirichlet data there. The singular 
integral equation for the Neumann data at any single point or a local patch on 
the boundary surface involves potential solution on a local hemisphere, which 
can be readily obtained with Feynman-Kac formula with the help of WOS 
sampling of the Brownian pathes. Numerical results validate the efficiency 
and accuracy of this method. 

The local BIE-WOS method of finding the DtN or NtD mapping can give 
a parallel algorithm for the solution of the Poisson equation with Dirichlet or 
Neumann boundary conditions. Firstly, we partition the whole boundary 
into a union of overlapping pathes Si namely, 

dn = UiSi, 

then, the local BIE-WOS can be used to find the DtN or NtD mapping over 
each patch Si independently in parallel. In principle, the computation of 
BIE-WOS over each patch can be done over one computing processor without 
need for communications with others, thus a high parallel scalability can be 
achieved. Secondly, the solution to the Poisson equation in the whole space 
can be found with the integral representation of (II. 3p with the help of one 
application of EMM [13] . 

There are several important research issues to be addressed before the 
BIE-WOS method can be used for large scale computation of Poisson or mod- 
ified Helmholtz equations, including (1) the NtD mapping problem, where the 
Neumann data is given on the boundary and the Dirichlet data is required. 
In this case, the Feynman-Kac formula derived in [18] can be used, which will 
involve reflecting Brownian pathes [20] with respect to the domain bound- 
ary. Efficient numerical implantation will have to be developed; (2) Modified 
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Helmholtz equation, even though the Feynman-Kac formula (|1.7p still applies, 
a survival factor will be introduced as the WOS samples the Brownian pathes 
and efficient way to use the Feynman-Kac formula will be addressed; and (3) 
the WOS scheme requires the computation of the distance between a Brown- 
ian particle and the boundary of the solution domain, efficient algorithms will 
be studied for the overall speed of the BIE-WOS method. 

The BIE-WOS based-parallel algorithm for solving Poisson or modified 
Helmholtz equation will have the following important features: 

• Non-iterative in construction and no need to solve any global linear 
system. 

• Stochastic in nature based on the fundamental link between Brownian 
motion and the solution of elliptic PDEs. 

• Massive parallelism suitable to large number of processors for large 
scale computing due to the random walk and local integral equation 
components of the algorithm. 

• No need for traditional finite element type surface or volume meshes. 

• Applicable to complex 3-D geometry with high accuracy treatment of 
domain boundary geometries. 

In comparison with traditional finite element and finite difference meth- 
ods, the BIE-WOS based solver is only suitable for Poisson and modified 
Helmholtz equations (due to the use of WOS-type sampling technique of the 
diffusion pathes) and its accuracy is limited to that of the Monte Carlo sam- 
pling technique, while the traditional grid based method can handle more 
general variable coefficients PDEs and achieve high accuracy. Nonetheless, as 
the Poisson and modified Helmholtz equations form the bulk computation of 
projection-type methods for incompressible flows and other important scien- 
tific computing, the progress in scalability of parallel BIE-WOS based-solvers 
will have large impact on the simulation capability of incompressible flow and 
engineering applications. 
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